{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0",
   "metadata": {
    "id": "0"
   },
   "source": [
    "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/segment-geospatial/blob/main/docs/workshops/IPPN_2024.ipynb)\n",
    "\n",
    "**Open Source Pipeline for UAS and satellite based High Throughput Phenotyping Applications - Part 1**\n",
    "\n",
    "This notebook is designed for workshop presented at the [International Plant Phenotyping Network (IPPN)](https://www.plant-phenotyping.org/ipps8) conference on October 7, 2024. Click the **Open in Colab** button above to run this notebook interactively in the cloud. For Part 2 of the workshop, please click [here](https://geemap.org/workshops/IPPN_2024/).\n",
    "\n",
    "-   Registration: <https://www.plant-phenotyping.org/index.php?index=935>\n",
    "-   Notebook: <https://samgeo.gishub.org/workshops/IPPN_2024>\n",
    "-   Earth Engine: <https://earthengine.google.com>\n",
    "-   Geemap: <https://geemap.org>\n",
    "-   Leafmap: <https://leafmap.org>\n",
    "-   Samgeo: <https://samgeo.gishub.org>\n",
    "-   Data to Science (D2S): <https://ps2.d2s.org>\n",
    "-   D2S Python API: <https://py.d2s.org>\n",
    "\n",
    "## Introduction\n",
    "\n",
    "Recent advances in sensor technology have revolutionized the assessment of crop health by providing fine spatial and high temporal resolutions at affordable costs. As plant scientists gain access to increasingly larger volumes of Unmanned Aerial Systems (UAS) and satellite High Throughput Phenotyping (HTP) data, there is a growing need to extract biologically informative and quantitative phenotypic information from the vast amount of freely available geospatial data. However, the lack of specialized software packages tailored for processing such data makes it challenging to develop transdisciplinary research collaboration around these data. This workshop aims to bridge the gap between big data and agricultural research scientists by providing training on an open-source online platform for managing big UAS HTP data known as Data to Science. Additionally, attendees will be introduced to powerful Python packages, namely leafmap and Leafmap, designed for the seamless integration and analysis of UAS and satellite images in various agricultural applications. By participating in this workshop, attendees will acquire the skills necessary to efficiently search, visualize, and analyze geospatial data within a Jupyter environment, even with minimal coding experience. The workshop provides a hands-on learning experience through practical examples and interactive exercises, enabling participants to enhance their proficiency and gain valuable insights into leveraging geospatial data for agricultural research purposes.\n",
    "\n",
    "## Agenda\n",
    "\n",
    "The main topics to be covered in this workshop include:\n",
    "\n",
    "* Create interactive maps using leafmap\n",
    "* Visualize drone imagery from D2S\n",
    "* Segment drone imagery using samgeo\n",
    "* Calculate zonal statistics from drone imagery\n",
    "* Visualize Earth Engine data\n",
    "* Create timelapse animations\n",
    "\n",
    "## Environment setup\n",
    "\n",
    "[![image](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/opengeos/segment-geospatial/blob/main/docs/workshops/IPPN_2024.ipynb)\n",
    "\n",
    "### Change Colab dark theme\n",
    "\n",
    "Currently, ipywidgets does not work well with Colab dark theme. Some of the leafmap widgets may not display properly in Colab dark theme.It is recommended that you change Colab to the light theme.\n",
    "\n",
    "![image](https://github.com/user-attachments/assets/dbb33f3c-084e-4913-9e06-2b66b8dcd4e3)\n",
    "\n",
    "\n",
    "### Change runtime type to GPU\n",
    "\n",
    "To speed up the processing, you can change the Colab runtime type to GPU. Go to the \"Runtime\" menu, select \"Change runtime type\", and choose \"T4 GPU\" from the \"Hardware accelerator\" dropdown menu.\n",
    "\n",
    "![image](https://github.com/user-attachments/assets/e92d2a19-0555-456d-b4be-36680c0af09f)\n",
    "\n",
    "\n",
    "### Install packages\n",
    "\n",
    "Uncomment the following code to install the required packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install -U \"leafmap[raster]\" segment-geospatial d2spy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %pip install numpy==1.26.4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3",
   "metadata": {
    "id": "3"
   },
   "source": [
    "### Import libraries\n",
    "\n",
    "Import the necessary libraries for this workshop."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4",
   "metadata": {
    "id": "4"
   },
   "outputs": [],
   "source": [
    "import leafmap"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5",
   "metadata": {
    "id": "9"
   },
   "source": [
    "## Creating interactive maps\n",
    "\n",
    "Let's create an interactive map using the `ipyleaflet` plotting backend. The [`leafmap.Map`](https://leafmap.org/leafmap/#leafmap.leafmap.m) class inherits the [`ipyleaflet.Map`](https://ipyleaflet.readthedocs.io/en/latest/map_and_basemaps/map.html) class. Therefore, you can use the same syntax to create an interactive map as you would with `ipyleaflet.Map`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7",
   "metadata": {
    "id": "11"
   },
   "source": [
    "To display it in a Jupyter notebook, simply ask for the object representation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8",
   "metadata": {},
   "outputs": [],
   "source": [
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9",
   "metadata": {
    "id": "13"
   },
   "source": [
    "To customize the map, you can specify various keyword arguments, such as `center` ([lat, lon]), `zoom`, `width`, and `height`. The default `width` is `100%`, which takes up the entire cell width of the Jupyter notebook. The `height` argument accepts a number or a string. If a number is provided, it represents the height of the map in pixels. If a string is provided, the string must be in the format of a number followed by `px`, e.g., `600px`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "10",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(center=[40, -100], zoom=4, height=\"600px\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "11",
   "metadata": {
    "id": "17"
   },
   "source": [
    "### Adding basemaps\n",
    "\n",
    "There are several ways to add basemaps to a map. You can specify the basemap to use in the `basemap` keyword argument when creating the map. Alternatively, you can add basemap layers to the map using the `add_basemap` method. leafmap has hundreds of built-in basemaps available that can be easily added to the map with only one line of code.\n",
    "\n",
    "Create a map by specifying the basemap to use as follows. For example, the `Esri.WorldImagery` basemap represents the Esri world imagery basemap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map(basemap=\"Esri.WorldImagery\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13",
   "metadata": {
    "id": "19"
   },
   "source": [
    "You can add as many basemaps as you like to the map. For example, the following code adds the `OpenTopoMap` basemap to the map above:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "14",
   "metadata": {
    "id": "20",
    "outputId": "b1fea60d-6a8c-4c9e-fc13-0e3525c4a917"
   },
   "outputs": [],
   "source": [
    "m.add_basemap(\"OpenTopoMap\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15",
   "metadata": {
    "id": "21"
   },
   "source": [
    "You can also add an XYZ tile layer to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16",
   "metadata": {
    "id": "22",
    "outputId": "c4f2108b-8a2b-4903-db1a-bcc0690a3651"
   },
   "outputs": [],
   "source": [
    "basemap_url = \"https://mt1.google.com/vt/lyrs=y&x={x}&y={y}&z={z}\"\n",
    "m.add_tile_layer(basemap_url, name=\"Hybrid\", attribution=\"Google\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17",
   "metadata": {
    "id": "23"
   },
   "source": [
    "You can also change basemaps interactively using the basemap GUI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_basemap_gui()\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19",
   "metadata": {
    "id": "25"
   },
   "source": [
    "## Visualizing Drone Imagery from D2S\n",
    "\n",
    "The Data to Science (D2S) platform (https://ps2.d2s.org) hosts a large collection of drone imagery that can be accessed through the D2S API (https://py.d2s.org). To visualize drone imagery from D2S, you need to [sign up](https://ps2.d2s.org/auth/register) for a free account on the D2S platform and obtain an API key."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20",
   "metadata": {
    "id": "de1a9ee4-d27e-49cb-9dcf-db32aab48ccc"
   },
   "source": [
    "### Login to D2S\n",
    "Login and connect to your D2S workspace in one go using the d2spy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21",
   "metadata": {},
   "outputs": [],
   "source": [
    "from d2spy.workspace import Workspace\n",
    "\n",
    "# Replace with URL to a D2S instance\n",
    "d2s_url = \"https://ps2.d2s.org\"\n",
    "\n",
    "# Login and connect to workspace with your email address\n",
    "workspace = Workspace.connect(d2s_url, \"workshop@d2s.org\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check for API key\n",
    "api_key = workspace.api_key\n",
    "if not api_key:\n",
    "    print(\n",
    "        \"No API key. Please request one from the D2S profile page and re-run this cell.\"\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from datetime import date\n",
    "\n",
    "os.environ[\"D2S_API_KEY\"] = api_key\n",
    "os.environ[\"TITILER_ENDPOINT\"] = \"https://tt.d2s.org\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24",
   "metadata": {
    "id": "TOx4ytjZKLOC"
   },
   "source": [
    "### Choose a project to work with\n",
    "\n",
    "The Workspace `get_projects` method will retrieve a collection of the projects your account can currently access on the D2S instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get list of all your projects\n",
    "projects = workspace.get_projects()\n",
    "for project in projects:\n",
    "    print(project)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "26",
   "metadata": {
    "id": "zTQJkyxaKxKn"
   },
   "source": [
    "The `projects` variable is a `ProjectCollection`. The collection can be filtered by either the project descriptions or titles using the methods `filter_by_title` or `filter_by_name`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "27",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of creating new collection of only projects with the keyword \"Citrus Orchard\" in the title\n",
    "filtered_projects = projects.filter_by_title(\"Citrus Orchard\")\n",
    "print(filtered_projects)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28",
   "metadata": {
    "id": "O6fdeffmLIN5"
   },
   "source": [
    "Now you can choose a specific project to work with. In this case, the filtered projects returned only one project, so we will use that project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "29",
   "metadata": {},
   "outputs": [],
   "source": [
    "project = filtered_projects[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30",
   "metadata": {},
   "source": [
    "### Get the project boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "31",
   "metadata": {
    "id": "yqo000pqLYn4"
   },
   "source": [
    "`get_project_boundary` method of the `Project` class will retrieve a GeoJSON object of the project boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get project boundary as Python dictionary in GeoJSON structure\n",
    "project_boundary = project.get_project_boundary()\n",
    "project_boundary"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33",
   "metadata": {
    "id": "OvFhMBNtLv9X"
   },
   "source": [
    "### Get project flights\n",
    "\n",
    "The `Project` `get_flights` method will retrieve a list of flights associated with the project."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get list of all flights for a project\n",
    "flights = project.get_flights()\n",
    "# Print first flight object (if one exists)\n",
    "for flight in flights:\n",
    "    print(flight)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35",
   "metadata": {
    "id": "eT8rrteTMM4z"
   },
   "source": [
    "### Filter flights by date\n",
    "\n",
    "The `flights` variable is a `FlightCollection`. The collection can be filtered by the acquisition date using the method `filter_by_date`. This method will return all flights with an acquisition date between the provided start and end dates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "36",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of creating new collection of only flights from June 2022\n",
    "filtered_flights = flights.filter_by_date(\n",
    "    start_date=date(2022, 6, 1), end_date=date(2022, 7, 1)\n",
    ")\n",
    "for flight in filtered_flights:\n",
    "    print(flight)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37",
   "metadata": {
    "id": "dOlbqBhRNX4e"
   },
   "source": [
    "Now, we can choose a flight from the filtered flight. Let's choose the flight on June 9, 2022."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38",
   "metadata": {},
   "outputs": [],
   "source": [
    "flight = filtered_flights[0]\n",
    "flight"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39",
   "metadata": {
    "id": "ngIYB3tANljF"
   },
   "source": [
    "### Get data products\n",
    "\n",
    "The Flight `get_data_products` method will retrieve a list of data products associated with the flight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get list of data products from a flight\n",
    "data_products = flight.get_data_products()\n",
    "\n",
    "for data_product in data_products:\n",
    "    print(data_product)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41",
   "metadata": {
    "id": "yiKBMsv7N4cy"
   },
   "source": [
    "The `data_products` variable is a `DataProductCollection`. The collection can be filtered by data type using the method `filter_by_data_type`. This method will return all data products that match the requested data type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of creating new collection of data products with the \"ortho\" data type\n",
    "ortho_data_products = data_products.filter_by_data_type(\"ortho\")\n",
    "print(ortho_data_products)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43",
   "metadata": {
    "id": "ZryOKEuROOnu"
   },
   "source": [
    "### Visualize ortho imagery\n",
    "\n",
    "Now we can grab the ortho URL to display it using leafmap."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_basemap(\"HYBRID\", show=False)\n",
    "ortho_data = ortho_data_products[0]\n",
    "ortho_url_202206 = ortho_data.url\n",
    "ortho_url_202206 = leafmap.d2s_tile(ortho_url_202206)\n",
    "m.add_cog_layer(ortho_url_202206, name=\"Ortho Imagery 202206\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45",
   "metadata": {
    "id": "33"
   },
   "source": [
    "### Visualize DSM\n",
    "\n",
    "Similarly, you can visualize the Digital Surface Model (DSM) from D2S using the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of creating new collection of data products with the \"dsm\" data type\n",
    "dsm_data_products = data_products.filter_by_data_type(\"dsm\")\n",
    "print(dsm_data_products)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47",
   "metadata": {},
   "outputs": [],
   "source": [
    "dsm_data = dsm_data_products[0]\n",
    "dsm_url_202206 = dsm_data.url\n",
    "dsm_url_202206 = leafmap.d2s_tile(dsm_url_202206)\n",
    "m.add_cog_layer(dsm_url_202206, colormap_name=\"terrain\", name=\"DSM 202206\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "48",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.cog_stats(dsm_url_202206)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49",
   "metadata": {
    "id": "35"
   },
   "source": [
    "Add a colorbar to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_colormap(cmap=\"terrain\", vmin=3, vmax=33, label=\"Elevation (m)\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51",
   "metadata": {},
   "source": [
    "### Visualize CHM\n",
    "\n",
    "Similarly, you can visualize the Canopy Height Model (CHM) from D2S using the code below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example of creating new collection of data products with the \"chm\" data type\n",
    "chm_data_products = data_products.filter_by_data_type(\"chm\")\n",
    "print(chm_data_products)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53",
   "metadata": {},
   "outputs": [],
   "source": [
    "chm_data = chm_data_products[0]\n",
    "chm_url_202206 = chm_data.url\n",
    "chm_url_202206 = leafmap.d2s_tile(chm_url_202206)\n",
    "m.add_cog_layer(chm_url_202206, colormap_name=\"jet\", name=\"CHM 202206\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "54",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.cog_stats(chm_url_202206)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_colormap(cmap=\"jet\", vmin=0, vmax=13, label=\"Elevation (m)\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56",
   "metadata": {
    "id": "37"
   },
   "source": [
    "Add the project boundary to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_geojson(project_boundary, layer_name=\"Project Boundary\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58",
   "metadata": {},
   "source": [
    "Add tree boundaries to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59",
   "metadata": {},
   "outputs": [],
   "source": [
    "map_layers = project.get_map_layers()\n",
    "tree_boundaries = map_layers[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_geojson(tree_boundaries, layer_name=\"Tree Boundaries\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61",
   "metadata": {},
   "source": [
    "### Get another flight\n",
    "\n",
    "Retrieve the Ortho data product for the December 2022 flight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "62",
   "metadata": {},
   "outputs": [],
   "source": [
    "filtered_flights = flights.filter_by_date(\n",
    "    start_date=date(2022, 12, 1), end_date=date(2022, 12, 31)\n",
    ")\n",
    "for flight in filtered_flights:\n",
    "    print(flight)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63",
   "metadata": {},
   "outputs": [],
   "source": [
    "flight_202212 = filtered_flights[0]\n",
    "data_products = flight_202212.get_data_products()\n",
    "ortho_data_products = data_products.filter_by_data_type(\"ortho\")\n",
    "ortho_data = ortho_data_products[0]\n",
    "ortho_url_202212 = ortho_data.url\n",
    "ortho_url_202212 = leafmap.d2s_tile(ortho_url_202212)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64",
   "metadata": {},
   "source": [
    "### Compare two ortho images\n",
    "\n",
    "Create a split map for comparing the 2022 and 2024 ortho images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipyleaflet import TileLayer\n",
    "\n",
    "m = leafmap.Map()\n",
    "left_layer = TileLayer(\n",
    "    url=leafmap.cog_tile(ortho_url_202206), max_zoom=30, name=\"2022-06 Ortho\"\n",
    ")\n",
    "right_layer = TileLayer(\n",
    "    url=leafmap.cog_tile(ortho_url_202212), max_zoom=30, name=\"2022-12 Ortho\"\n",
    ")\n",
    "m.split_map(left_layer, right_layer, left_label=\"2022-06\", right_label=\"2022-12\")\n",
    "m.set_center(-97.955281, 26.165595, 18)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "66",
   "metadata": {},
   "source": [
    "## Download data from D2S\n",
    "\n",
    "Read the ortho image from D2S as a DataArray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67",
   "metadata": {},
   "outputs": [],
   "source": [
    "import rioxarray as rxr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = rxr.open_rasterio(ortho_url_202206)\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "69",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_cog_layer(ortho_url_202206, name=\"Ortho Imagery 202206\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70",
   "metadata": {},
   "source": [
    "Draw an area of interest (AOI) on the map. If an AOI is not provided, a default AOI will be used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71",
   "metadata": {},
   "outputs": [],
   "source": [
    "if m.user_roi is not None:\n",
    "    bbox = m.user_roi_bounds()\n",
    "else:\n",
    "    bbox = [-97.956252, 26.165315, -97.954992, 26.165883]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72",
   "metadata": {},
   "outputs": [],
   "source": [
    "geojson = leafmap.bbox_to_geojson(bbox)\n",
    "gdf = leafmap.geojson_to_gdf(geojson)\n",
    "m.add_gdf(gdf, layer_name=\"AOI\", info_mode=None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "73",
   "metadata": {},
   "outputs": [],
   "source": [
    "crs = data.rio.crs.to_string()\n",
    "print(crs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = gdf.to_crs(crs)\n",
    "print(gdf.crs)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75",
   "metadata": {},
   "source": [
    "Resample the ortho imagery from 1 cm to 10 cm resolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "76",
   "metadata": {},
   "outputs": [],
   "source": [
    "resampled_data = data.rio.reproject(crs, resolution=(0.1, 0.1))\n",
    "resampled_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77",
   "metadata": {},
   "source": [
    "Clip the ortho image to the AOI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipped_data = resampled_data.rio.clip(gdf.geometry, gdf.crs)\n",
    "clipped_data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79",
   "metadata": {},
   "source": [
    "Save the clipped ortho image to a GeoTIFF file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80",
   "metadata": {},
   "outputs": [],
   "source": [
    "image = \"ortho_image_202206.tif\"\n",
    "clipped_data.sel(band=[1, 2, 3]).rio.to_raster(image)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81",
   "metadata": {},
   "source": [
    "Read the CHM dataset from D2S as a DataArray."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82",
   "metadata": {},
   "outputs": [],
   "source": [
    "chm_data = rxr.open_rasterio(chm_url_202206)\n",
    "chm_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83",
   "metadata": {},
   "outputs": [],
   "source": [
    "resampled_chm_data = chm_data.rio.reproject_match(resampled_data)\n",
    "resampled_chm_data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84",
   "metadata": {},
   "outputs": [],
   "source": [
    "clipped_chm_data = resampled_chm_data.rio.clip(gdf.geometry, gdf.crs)\n",
    "clipped_chm_data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85",
   "metadata": {},
   "outputs": [],
   "source": [
    "chm_image = \"chm_202206.tif\"\n",
    "clipped_chm_data.sel(band=[1]).rio.to_raster(chm_image)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86",
   "metadata": {},
   "source": [
    "Visualize the clipped ortho image."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_raster(image, layer_name=\"Ortho Image 202206\")\n",
    "m.add_geojson(tree_boundaries, layer_name=\"Tree Boundaries\")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "88",
   "metadata": {},
   "source": [
    "## Segmenting Drone Imagery using Samgeo\n",
    "\n",
    "### Initialize SAM class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89",
   "metadata": {},
   "outputs": [],
   "source": [
    "from samgeo import SamGeo, SamGeo2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "90",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2 = SamGeo2(model_id=\"sam2-hiera-large\", automatic=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91",
   "metadata": {},
   "source": [
    "### Automatic mask generation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "92",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.generate(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.save_masks(output=\"masks.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.show_masks(cmap=\"binary_r\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95",
   "metadata": {},
   "source": [
    "![image](https://github.com/user-attachments/assets/fec9ad33-338e-4d1f-95e5-5f7b58e15edb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.show_masks(cmap=\"jet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "97",
   "metadata": {},
   "source": [
    "![image](https://github.com/user-attachments/assets/32036cec-96a2-46c7-b948-a859925d8600)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98",
   "metadata": {},
   "source": [
    "Show the object annotations (objects with random color) on the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.show_anns(axis=\"off\", alpha=0.7, output=\"annotations.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "100",
   "metadata": {},
   "source": [
    "![image](https://github.com/user-attachments/assets/8a05f797-5a08-4bb2-b6b3-ee8f82d8bd32)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "101",
   "metadata": {},
   "source": [
    "Compare images with a slider."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "102",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.image_comparison(\n",
    "    image,\n",
    "    \"annotations.tif\",\n",
    "    label1=\"Drone Imagery\",\n",
    "    label2=\"Image Segmentation\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "103",
   "metadata": {},
   "source": [
    "Add segmentation result to the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "104",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\"masks.tif\", colormap=\"jet\", layer_name=\"Masks\", nodata=0, opacity=0.7)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "105",
   "metadata": {},
   "source": [
    "Convert the object masks to vector format, such as GeoPackage, Shapefile, or GeoJSON."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "106",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.raster_to_vector(\"masks.tif\", \"masks.gpkg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "107",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_vector(\"masks.gpkg\", layer_name=\"Objects\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "108",
   "metadata": {},
   "source": [
    "### Automatic mask generation options\n",
    "\n",
    "There are several tunable parameters in automatic mask generation that control how densely points are sampled and what the thresholds are for removing low quality or duplicate masks. Additionally, generation can be automatically run on crops of the image to get improved performance on smaller objects, and post-processing can remove stray pixels and holes. Here is an example configuration that samples more masks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "109",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2 = SamGeo2(\n",
    "    model_id=\"sam2-hiera-large\",\n",
    "    apply_postprocessing=False,\n",
    "    points_per_side=64,\n",
    "    points_per_batch=128,\n",
    "    pred_iou_thresh=0.7,\n",
    "    stability_score_thresh=0.92,\n",
    "    stability_score_offset=0.7,\n",
    "    crop_n_layers=1,\n",
    "    box_nms_thresh=0.7,\n",
    "    crop_n_points_downscale_factor=2,\n",
    "    min_mask_region_area=25,\n",
    "    use_m2m=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "110",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.generate(image, output=\"masks2.tif\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "111",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.show_masks(cmap=\"jet\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "112",
   "metadata": {},
   "source": [
    "![image](https://github.com/user-attachments/assets/ff07bff9-f0d1-41c6-a641-ce9f4734fd61)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "113",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam2.show_anns(axis=\"off\", alpha=0.7, output=\"annotations2.tif\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "114",
   "metadata": {},
   "source": [
    "![image](https://github.com/user-attachments/assets/1e1df62c-ec44-4c5e-aa5a-4025d078ff6b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "115",
   "metadata": {},
   "outputs": [],
   "source": [
    "leafmap.image_comparison(\n",
    "    image,\n",
    "    \"annotations2.tif\",\n",
    "    label1=\"Image\",\n",
    "    label2=\"Image Segmentation\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "116",
   "metadata": {},
   "source": [
    "Remove small objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "117",
   "metadata": {},
   "outputs": [],
   "source": [
    "da, gdf = sam2.region_groups(\n",
    "    \"masks2.tif\",\n",
    "    connectivity=1,\n",
    "    min_size=10,\n",
    "    max_size=2000,\n",
    "    intensity_image=\"chm_202206.tif\",\n",
    "    out_image=\"objects.tif\",\n",
    "    out_csv=\"objects.csv\",\n",
    "    out_vector=\"objects.gpkg\",\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "118",
   "metadata": {},
   "source": [
    "### Using box prompts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "119",
   "metadata": {},
   "outputs": [],
   "source": [
    "gdf = leafmap.geojson_to_gdf(tree_boundaries)\n",
    "gdf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "120",
   "metadata": {},
   "outputs": [],
   "source": [
    "geojson = \"tree_boundaries.geojson\"\n",
    "gdf.to_file(geojson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "121",
   "metadata": {},
   "outputs": [],
   "source": [
    "m = leafmap.Map()\n",
    "m.add_raster(image, layer_name=\"image\")\n",
    "style = {\n",
    "    \"color\": \"#ffff00\",\n",
    "    \"weight\": 2,\n",
    "    \"fillColor\": \"#7c4185\",\n",
    "    \"fillOpacity\": 0,\n",
    "}\n",
    "m.add_vector(geojson, style=style, zoom_to_layer=True, layer_name=\"Bounding boxes\")\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "122",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam = SamGeo(\n",
    "    model_type=\"vit_h\",\n",
    "    automatic=False,\n",
    "    sam_kwargs=None,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "123",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam.set_image(image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "124",
   "metadata": {},
   "outputs": [],
   "source": [
    "sam.predict(\n",
    "    boxes=geojson, point_crs=\"EPSG:4326\", output=\"tree_masks.tif\", dtype=\"uint16\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "125",
   "metadata": {},
   "outputs": [],
   "source": [
    "m.add_raster(\n",
    "    \"tree_masks.tif\", cmap=\"jet\", nodata=0, opacity=0.5, layer_name=\"Tree masks\"\n",
    ")\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "126",
   "metadata": {},
   "source": [
    "![image](https://github.com/user-attachments/assets/0a9b794d-cf94-4e5c-9134-8799e67d338c)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "sam",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
